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WEAK  SHOCK  WAVES  PROPAGATING  IN  FLUID-GAS 
MIXTURES 


1 .  INTRODUCTION 


Suspensions  of  solid,  liquid  or  gaseous  particles  in  an  ambient 
fluid  are  common  both  in  nature  and  industrial  processes.  Interestinq 
military  examples  are  fuel-air  explosives,  aerated  liquid  explosives,  foams 
(Dewey  1981,  1983)  and  bubble  curtains  (Mader  1982).  fliere  is  some  interest 
therefore  to  understand  the  basic  physical  mechanisms  that  occur  when  a 
liquid-gas  mixture  is  shock  loaded. 

It  is  well  known  that  the  presence  of  gas  bubbles  in  an  homogeneous 
explosive  significantly  increases  the  sensitivity  of  the  explosive,  Campbell 
et.al.  (1961),  Evans  et.  al.  (1962)  and  Mader  (1965).  Ttieir  results  show  that 
when  bubbles  in  the  explosive  are  heated  significantly  above  the  temperature 
of  the  liquid  due  to  impact  with  the  walls  of  the  bubble,  hot  spots  form  in 
the  explosive  and  detonation  can  result. 


The  attenuation  of  shock  waves  by  foams  is  of  particular  interest 
where  protection  from  blast  waves  is  needed.  Liquid  foams  have  physical 
characteristics  such  as  low  velocity  of  propagation  of  sound,  hiqh  diffusion 
losses,  and  high  heat  capacity  losses  making  them  suitable  as  blast 
suppressants. 


Shock  waves  in  fluids  containing  qas  bubbles  have  been  considered  to 
some  extent  by  Crespo  (1969)  and  Campbell  and  Pitcher  (1957)  and  the 
references  contained  therein.  Crespo  considers  the  propagation  of 
infinitesimal  sound  waves  in  a  fluid  containing  gas  bubbles  and  the  structure 
of  steady  shock  waves.  Campbell  and  Pitcher  have  derived  the  Huqoniot 
relations  (i.e.  relations  that  express  conservation  of  mass  and  momentum 
across  the  shock)  for  a  normal  steady  shock  in  a  fluid  containing  gas  bubbles, 
neglecting  temperature  rises  across  the  shock.  They  have  also  reported 
experimental  measurements  of  shock  waves  in  bubble-liouid  mixtures  using  a 
small  gas-liquid  shock  tube. 
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In  the  following  paper  a  simple  two-phase  flow  model  is  developed 
for  the  shock  transition  region  of  a  weak  shock  wave  propagating  into  a  fluid 
containing  spherical  gas  bubbles  whose  size  is  small  compared  to  the  shock 
transition  zone.  lhe  shock  wave  is  generated  by  the  impulsive  motion  of  a 
piston.  Assumptions  made  in  deriving  an  analytical  solution  are  similar  to 
those  of  Moran  and  Shen  (1966),  Campbell  and  Pitcher  (1977),  Crespo  (1969)  and 
Van  Wijngaarden  (1970).  They  are  (i)  the  fluid  is  incompressible  and  the  gas 
follows  the  perfect-gas  law,  (ii)  boundary  layer  effects  can  be  neglected, 
(iii)  the  only  viscous  effects  are  due  to  the  gas  bubble  draq  and  the  only 
heat  transfer  is  that  between  the  surrounding  liguid  and  qas  bubbles,  (iv)  the 
thermal  motion  of  the  gas  bubbles  is  negligible  and  so  does  not  contribute  to 
the  pressure  and  (v)  the  shock  wave  moves  with  the  velocity  of  the  piston. 

The  shock  structure  is  investigated  in  the  context  of  the  linear 
theory  i.e.,  the  near  field  solution.  The  methodology  is  to  expand  all 
dependent  state  variables  (which  can  be  collectively  designated  *)  in  the 
Navier-Stokes  equations  describing  the  flow,  as  f  +  e  +  ••••» 

and  equate  terms  of  order  0(c)  or  higher  to  zero.  The  procedure  to  be 
followed  in  solving  this  set  of  coupled  linear  partial  differential  equations 
with  assigned  boundary  and  initial  conditions  is  to  use  integral  transform 
methods.  Hie  solution  gives  an  approximation  to  the  near-field  fluid  flow. 


2.  BASIC  EQUATIONS 


The  shock  transition  for  a  weak  shock  wave  propagating  through  an 
inviscid  and  incompressible  fluid  containing  gas  bubbles  is  shown  in  figure  1 


u  «u0«u 

g  1  p 


Figure  1.  Shock  wave  propagating  through  a  fluid-gas  mixture. 


The  conditions  of  the  gas  are  described  by  its  velocity  (u  )  and  bv 
two  state  variables  density  (p  )  and  temperature  (T  ).  The  suffix  r  or  g  to 


a  symbol  indicates  that  it  refirs  to  the  liquid  or  gas  component.  Other 


state  variables  such  as  pressure  (p  )  or  the  speed  of  sound  may  be  introduced 
through  appropriate  thermodynamic  relations  such  as  an  equation  of  state  of 
the  gas.  Since  the  radius  of  the  gas  bubbles  are  assumed  small  compared  to  the 
thickness  of  the  shock,  the  pressure  of  the  gas  (P  )  can  be  assumed  equal  to 
that  of  the  liquid  (P#).  The  state  variables  of  the  liquid  are  characterised 


by  its  velocity  (u,),  temperature  (T. )  and  volume  fraction  of  qas  in  the 
fluid  (0). 
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In  a  cartesian  coordinate  system,  the  one-dimensional  tine  dependent 

equations  describing  the  conservation  of  mass,  momentum  and  energy  of  the 

fluid-gas  mixture  are  given  as  follows.  Let  the  density  of  the  mixture  p  be 

in 

given  in  terms  of  the  densities  of  the  two  components  by:  =  (1-B)p  +  Bp  . 

Assuming  no  mass  exchange  between  the  phases,  the  continuity  equations  q 
expressed  in  terms  of  the  volume  fraction  of  gas  B  and  volume  fraction  of 
liquid  (1-8)  are: 


3t(1-B)  +  3x(l-6)uJt  -  0, 

for  the  liquid  and 


(2.1 ) 


3  (Bp  )  +  3  (Bp  u  )  *  0, 
t  g  x  g  g 


(2.2) 


for  the  gas  bubbles,  where  the  density  of  the  liquid  in  (2.1)  has  been 
cancelled  out  because  the  fluid  is  assumed  incompressible. 

Hie  equation  of  momentum  of  the  fluid-gas  mixture  is  given  by 


where 


•a- 


—  :  ■*  3.  ♦  u, . .  3  . 

Dt(j)  t  (j)  x 


(2.3) 


is  the  differential  following  the  motion  of  the  fluid,  the  index  j  can  be 
either  g  or  l  for  gas  or  liquid  respectively. 


The  momentum  equation  for  a  volume  of  gas  can  be  written  following 
Crespo  (1969)  as 


3  P 
x  g 


+ 


X 
v  ' 


(2.4a) 


where  u^:  »  u^  -  ug,  f  is  the  internal  force,  V  is  the  volume  of  a  qas 
bubble,  and  T  a  function  which  depends  on  B.  For  spherical  bubbles 
and  8  small,  T  *  1 .  The  force  on  a  qas  bubble  movinq  throuqh  the  fluid  has 
been  considered  by  Van  Wijngaarden  (1970)  and  Levich  (1962).  For  hiqh 
Reynolds  numbers  the  relative  motion  of  the  gas  bubble  can  be  determined  from 
the  irrotational  inviscid  flow  around  the  bubble,  since  the  gas  bubble  is 
free  to  move  there  is  no  constraint  on  the  tangential  velocity  of  the  fluid  at 
the  surface  of  the  bubble.  Hence  no  velocity  boundary  layer  exists.  Thus  the 
frictional  force  acting  on  the  bubble  can  be  taken  to  approximate  a  Stokes  law 
i.e,,  the  drag  force  on  a  spherical  bubble  is  proportional  to  the  relative 
fluid  velocity: 
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6ifr  u  u.  , 

q  1  lq 


where  is  the  dynamic^ viscosity  of  the  fluid  and  rq  the  initial  radius  of  a 
gas  sphere.  The  tern  —  appearing  in  (2.4a)  can  now  be  expressed  as  follows: 


9 

2  UtqP 


4wp  2/3 


tV-3T! 


where  is  the  kineaatic  viscosity  of  the  fluid.  This  equation  represents 
the  dissipative  force  per  unit  volume  which  in  terns  of  the  draq  coefficient 
Cp  can  be  written  as 


where 


3  it  , 

a  8  r  uiglUlg'CD' 

9 


1 2v  4xp  1/3 


In  terns  of  the  monentm  relaxation  time  a  : 

« 

a  — E —  - 

.  ewr^  9  wt 


the  nomentum  equation  of  the  gas  can  be  expressed  as 


3  P 

x  g 


.  Du.  u 

1  „  lq  ig 

2  Pt  Dt,  ,  g  a 

(g)  m 


(2.4b) 


where  the  second  tern  of  (2.4b)  represents  a  non-dissipative  force  per  unit 
voluae  due  to  the  mass  of  the  bubble. 


The  equation  of  energy  of  the  gas  can  be  written  as 


DP  Y  P  Dp 

—5  _  -2-3.  —3  m 

Dt,  I  P  Dt 
(g)  g  (g) 


(Yq  -  Do, 


(7.5) 


where  Q  is  the  heat  transfer  rate  and  y  the  ratio  of  specific  heats  of  the 
gas.  Crespo  (1969)  gives  p  as  q 
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- *  r  o.(T  -  T  )Nu, 

m  g  Jt  o  g 


where  <J  and  TQ  are  the  thermal  conductivity  and  the  initial  temperature  of 
the  liquid  respectively.  Nu  is  the  Nusselt  number.  For  weak  wave  analysis 

we  assume  there  is  only  pure  heat  conduction,  thus  Nu  -  2  and 


4s p  2/3 

3(^  V*. 


.  T  )  . 

to  g 


Hie  final  equation  is  the  simple  equation  of  state  of  the  gas: 


P  »  p  RT  .  (2.6) 

9  9  9 


Equations  (2.1)  through  to  (2.6)  give  six  simultaneous  eguations  in  x  and  t 
for  the  fluid  flow  variables  ♦  »  (P^,  u^,  ug,  T^,  p^,  $).  These  must  be 
solved  for  the  spatial  and  temporal  dependences.  Where  a  shock  wave 
propagates  down  a  semi-infinite  shock  tube,  the  initial  conditions  are 


u,  *  u  =0;  x  >  0, 
*  9 


P  . 

o  g 


po'  Ti 


V 


and  the  boundary  conditions  are 


U.  ■  U  ,  t  >  0,  x«ut. 
Ip  p 


u 

9 


0,  u.  ♦  0,  P,  P_ 


as  x  +  ® 


It  is  necessary  to  introduce  a  dimensionless  perturbation 
parameter  e  which  is  a  measure  of  the  shock  strength. 


e 


where  U  i3  the  piston  velocity  and  aQ  is  the  equilibrium  speed  of  sound  in 
the  undisturbed  flow.  All  the  dependent  state  variables  \J<  can  be  expanded  in 
terms  of  t  as  (Van  Dyke  1975) : 


(0) 


4- 


(D 


2,  (2) 

e  + - , 


and  in  particular 
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* - ) 


(2.7a) 


ut  ■  “o1'!"  *  “I21  ♦  eM31 


»  -  «  (u">  .  eu12’  *  .V3’.  -  -  -I, 

g  o  g  g  g 


(2.7b) 


0  -  0  (1  +  £0(1)  +  e20(2)+  -  -  -), 
o 


(2.7c) 


p  «  P  (1  +  ep(1)  +  e2p(2)+  -  -  -), 
g  o  g  g  ' 


(2. 7d) 


T  -  T  (1  +  6T(1 )  +  £2T(2)+ 
go  g  g 


(2.7e) 


p  -  p  (1  +  ep*1  *  +  C2p(2,+ - ), 

go  g  g 


( 2. 7f ) 


where  represents  the  state  variable  in  the  Initial  stationary 
state  and  ^in+1 ) /^(n)  is  0£  order  0(e). 

The  perturbation  scheme  is  illustrated  with  the  continuity  equation 
of  the  liquid  (2.1).  After  a  change  of  variable  from  (x,t)  to  the 
dimensionless  variables  (£,x): 


5 


the  continuity  equation  (2.1)  becomes 


T 


t. 


a  3(1 -0)  +  3,(1 -0)u ,  -  0. 

O  T  <;  t 


(2.8) 


Upon  substitution  of  the  perturbation  expansions  (2.7)  in  equation  (2.8)  we 
obtain  to  order  0(e) s 


0  3  0(1)  +  (1-0  )3ru[n  -  0. 
ox  o  £  t 


(2.9a) 


Proceeding  in  a  like  manner  for  equations  (2.2)  to  (2.6)  we  obtain  to  order 
0(e)  a  system  of  partial  differential  equations  where  all  dependent 
perturbation  variables  are  of  first  order: 
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■iafc£ 


■ 

VV 


V*  ■  * 


£ 


S  tn(1) 

Vpg 


8(,)) 


v(,) 

5  g 


o, 


(2.9b) 


3  P(1  )  -  Y  3  P(1  )  +  =  0, 

T  g  g  T  g  g 


(2.9c) 


3 ,P(1)  +  T  3_(u(1)  -«<’>)  ♦  (u(1)  -  u‘1})a*  =  0, 

<;  g  JTg  *  g  t 


( 2. 9d ) 


(1-8  lau''1  +  3fpn)  -  0, 
o  t  l  C  g 


(2.9e) 


_(1)  M)  _(1)  . 

p  -  p  -  T  »  0, 

g  g  g 


( 2. 9f ) 


where  a*  *  9v^(2UpTg)-1  is  the  Reynolds  number  of  a  typical  gas 

bubble  and  a  -  Ygrg(UpCkp) .  Here  is  the  thermal  relaxation  time  of  a  gas 

bubble  defined  by  otj.  »  mcp(4*rga£)-1. 

Hie  system  of  equations  (2.9)  can  be  solved  by  taking  the  Laplace 
transform  of  *  with  respect  to  the  time  variable  T: 


*(5,3) 


f  e  3T*(C,T)dT, 

o 


(2.10) 


Hie  operation  of  multiplying  equation  (2.9)  by  the  kernel  of  the  Laplace 
transform  and  integration  with  respect  to  t  over  the  interval  (0,*»)  together 
with  (2.10)  and  the  initial  conditions,  give  set  of  six  coupled  auxiliary 
equations.  Decoupling  this  set  of  equations  for  the  state  variables  *  leads 
to  the  algebraic  equation: 


-(1)  2  -(1) 
Ut,C5  '  3  h(S,Ul 


0, 


(2.11) 


where  the  coefficient  h(s)  is: 


(a  +  s) ( 1  -  6  ) ( Ts/2  +  a*J0 
.  ,  .  _ o _ o _ 

'  “  (T  s  +  a) (s( T/2  +  0  (1-0  ))  +  a*)’ 

g  o  o 


This  procedure  has  reduced  the  problem  to  the  solution  of  an  ordinary 
differential  equation  (2.11),  provided  s  *  -ay  -1  and 

s  *  -a*(r/2  ♦  0O(1-0O))“1.  The  solution  of  (2.11),  for  5  >  0  satisfying  the 
boundary  conditions  is 
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(2.1 2a) 


-( 1 )  1  -s/h(s)£  , ,  . 

u  (C,s)  =  —  e  .  (2.12a) 

i  s 

The  inversion  to  uj1^(5,t)  is  now  accomplished  by  use  of  the  inversion  formula 
for  Laplace  transforms  and  gives 


(1),,  . 

(C,T) 


-V  / 

2111  z 


exp ( st  -  s/h(s  11) 


(2.1 2b) 


where  Z  is  the  path  of  integration,  such  that  Re(s)  is  constant  and  is 
situated  to  the  right  of  all  singularities.  An  analytic  solution  to  equation 
(2.12b)  can  be  derived  for  certain  limiting  and  important  cases.  For  s  larqe, 
we  have  the  case  corresponding  to  high  frequencies  and  consequently  this 
approximation  is  valid  when  the  high  frequency  waves  dominate  ie. ,  t  small,  or 
near  discontinuities  in  the  wave  form.  For  s  small,  corresponding  to  T  large, 
the  high  frequency  waves  are  attenuated  so  that  lower  frequency  waves 
dominate.  For  s  small  the  expansion  for  h(s)  is 


h(s)  =  /5  -sa/b  +  0(s2), 


where  the  dimensionless  constants  A  and  B  are  defined  by 


A-U-2 

2  v  a 


y  -  i  -2  i, 

-3 —  +  -2--),  b  -  (e  d-6  ))"  '2 . 

a  a*  ‘  v  o  o  ' 


The  constant  B  represents  a  dimensionless  velocity  at  which  the  wave 
propagates  and  has  a  minimum  value  when  8  ”  ^2  •  T*le  constant  A  incorporates 

the  effects  of  viscosity  and  thermal  conductivity. 


3.  SOLUTION  FOR  LARGE  TIMES 


Using  the  method  of  steepest  descent  the  above  integral  can  be 
evaluated  (Van  Dyke  (1975))  and  the  solution  in  non-dimensional  form,  is  qiven 
by 

(D  1  ,5-RT 

u  { 1 , T  )  «  —  erfc( - — )  .  (3.1) 

1  2  2B/AT 


When  the  above  solution  i3  substituted  into  (2.9)  we  find 


(1)  H)  (1) 

u.  =  u  ,  T  »  0, 

l  q  g 


and  for  the  pressure  and  density: 


i  ;  \m>  IV  UJjllUlhMi" 


g 


v'  . 


g 


(3.3) 


It  can  be  seen  from  (3.1)  that  this  wave  travels  with 
velocity  BaQ  in  an  (x,t)  coordinate  system,  and  diffuses  due  to  the  combined 
effects  of  viscosity  and  thermal  conductivity  with  a  characteristic  diffusion 
width  given  by: 


5  -  BT  =  2B/AT, 


°r 

r  At  V2 

x  -  a  Bt  =  2a  b(— "J — ) 
o  o  v  U  ' 

P 


The  picture  of  the  fluid  and  gas  emerging  from  the  above  analysis  is  that, 
provided  t  is  large,  the  bubbles  behave  isothermally  (since  T  =  Tq  from 
(2.7e))  and  propagate  with  the  velocity  of  the  fluid  (u„  =  u  ).  The  wave 

X  Q 

described  by  (3.1)  has  a  further  property  that  as  A  +  0,  the  wave  becomes 
steeper  and  when  it  becomes  discontinuous  then  the  solution  represents  a  wave 
for  which  no  dissipative  processes  are  present. 


I 

4.  SOLUTION  FOR  SMALL  TIMES 


An  approximate  evaluation  of  the  integral  (2.12b)  corresponding  to 
the  hiqh  frequency  waves  (i.e.,  s  large,  small  time)  can  be  found  by  expanding 
s  large  in  h(s)  to  give: 


-1  2 
h(s)  =  a  +  abs  +  0(1 /s  ) 


where 


=  (' 


Y  (Tb  +  2) 

q 


~)  ^  and  b  *  (- 


a  (Y  -  D 


4aY  a* 


r)- 


Substitution  of  h(s)  into  (2.12a)  and  integration  leads  to  a  liquid  dynamic 
wave  whose  velocity  is  given  by 


(1)  -abC 

Ujj  (C  )  -  H(x  -  £a)e  , 


(4.1  ) 


an«i  moreover 
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u*(5,t)  -  ea  H(t  -  £a)e  +  0(e  )  , 

where  H(x)  is  the  Heaviside  step  function  defined  as  H(x)  =  1  for  x>0  and 
H(x)  »  0  for  x<0.  the  above  equation  indicates  a  wave  which  travels  at  a 
speed  a,  and  decays  exponentially  with  a  decay  length  x  =  r  aQ(abOp)_1 , 
where  t  is  the  mean  free  path  for  transfer  of  momentum  by  the  gas ^bubbles . 
the  first  order  perturbations  in  pressure,  density,  temperature  and  velocity 
of  the  gas  are  respectively: 


P(1)  -  (1-0  i  )a"1H(T-?a>e“*b5  -  Y^V15  ,  (4.2) 

go  g  g 

T(1)  -  (Y  -1)0-0  )Y“1H(T-C»)«"*b5,  (4.3) 

g  g  o  g 

u<15  -  2(0 -0Q)  +  V2  1*) Th(  T-Ca)e~ab^  . 


Qualitatively,  the  gas  bubbles  move  faster  than  the  fluid  since 
u^}<  Ug1 '  i.e.,  uA<u  .  Prom  (4.3)  and  (2.7e)  we  find  that  the  temperature  of 
the  gasyis  greater  than  the  fluid  (i.e.,  T  >Te).  the  wave  for  small  time, 
propagates  at  the  speed  a,  and  decays  exponentially  due  to  the  effects  of 
viscosity  and  thermal  conductivity.  The  process  occurs  i sen tropically . 


5.  DISCUSSION 


A  theoretical  model  has  been  developed  describing  a  two-phase  liquid 
gas  mixture  in  one  dimension.  In  this  model,  the  Navier-Stokes  equations  are 
firstly  linearized  then  decoupled  and  finally  solved  by  Laplace  transform 
technique.  He  find  that  for  small  times  high  frequency  waves  dominate  and  the 
gas  bubbles  in  the  liquid  propagate  faster  than  the  fluid.  Thermodynamically 
the  bubbles  behave  isentropically . 


In  the  case  where  times  are  large,  it  is  found  that  the  high 
frequency  waves  are  attenuated  so  that  the  lower  frequency  waves  dominate. 
The  gas-bubbles  in  this  case  behave  isotherraally  with  the  shock  wave 
approaching  a  steady  state.  For  the  low  frequency  case,  the  width  of  the 
transition  zo.-.e  or  shock  thickness  5,  based  on  the  maximum  slope: 

U 


(‘3T) 
at,  max 

can  be  calculated  and  leads  to  6  increasing  indefinitely  (i.e.,  as  t  ♦  ") 
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as  /"t.  wpor  a  steady  weak  shock  wave  it  is  well  known  that  6  -  0(e_1),  so  as 
t  ♦  •,  we  would  expect  the  linearized  solution  to  break  down  when 
t  =»  0 (e"2).  That  is  when  t  is  large  compared  to  e~2,  the  linearized  solution 
yields  an  overthick  shock  wave.  For  this  case  a  uniformly  valid  solution  can 
be  found  using  singular  perturbation  methods,  and  will  be  discussed  in  a 
following  paper. 
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